arXiv:astro-ph/9808354vl 31 Aug 1998 


Gravitational Leasing by Power-Law Mass Distributions: A Fast 

and Exact Series Approach 


Kyu-Hyun ChaeQ, Valery K. Khersonsky^, and David A. Turnshek^ 


Received 


accepted 


^Department of Physics & Astronomy, University of Pittsburgh, Pittsburgh, PA 15260 





ABSTRACT 


We present an analytical formulation of gravitational leasing using familiar triaxial 
power-law mass distributions, where the 3-dimensional mass density is given by 
p{X,Y,Z) = po [l + + (t)^ deflection angle and magnification 

factor are obtained analytically as Fourier series. We give the exact expressions for the deflection 
angle and magnification factor. The formulae for the deflection angle and magnification factor 
given in this paper will be useful for numerical studies of observed lens systems. An application 
of our results to the Einstein Cross can be found in Chae, Turnshek, & Khersonsky (1998). 
Our series approach can be viewed as a user-friendly and efficient method to calculate lensing 
properties that is better than the more conventional approaches, e.g., numerical integrations, 
multipole expansions. 


Subject headings: cosmology: gravitational lensing — galaxies: structure — methods: analytical 



1. Introduction 


Gravitational lenses in nature usually consist of a primary lensing galaxy and, in many cases, 
perturbation field caused by, for example, nearby galaxies and clusters of galaxies along the sight-line 
(see Keeton, Kochanek, & Seljak 1997). Thus, modeling an observed lens requires the construction of a 
realistic mass model, including accounting for any perturbation that may be present. This paper developes 
a (semi-)analytical method to calculate lensing due to the primary mass distribution. 

Lensing objects (e.g. galaxies or galaxy clusters) show, in general, 3-dimensional mass distributions 
which are non-spherically symmetric. Elliptical galaxies are known to be triaxial in shape, as well as 
oblate and prolate (see, e.g., Ryden 1992). Spiral galaxies can be considered flattened spheroids. Clusters 
of galaxies are also known to be non-spherical in shape. Gravitational lensing by these concentrations of 
mass in the universe is difficult to deal with mathematically. For a projected surface density of elliptical 
shape Schramm (1990) introduced a relatively simple formulation equivalent to the early formulation by 
Bourassa, Kantowski, & Norton (1973) and Bourassa & Kantowski (1975, see also Bray 1984). In the 
above formulations the deflections are given as integrals. Thus, one can deal with lensing by elliptical mass 
distributions through numerical integrations. While numerical integration is a way of dealing with lensing 
by elliptical mass distributions, it has not been preferred by most researchers. 

To “simulate” real lenses in nature while avoiding mathematical complexities or costly numerical 
integrations, several lens models have been considered in the literature. First, circularly symmetric lenses 
perturbed by an external shear term (also called a quadrupole term) were used (e.g. Chang & Refsdal 
1979, 1984; Kovner 1987a; Kochanek 1991; Wambsganss & Paczyhski 1994). In addition, Kochanek (1991) 
considered other types of perturbations (internal & mixed). This approach has the great advantage of 
mathematical simplicity. It could be a good approximation of a real lens in some situations. However, this 
approach is somewhat artificial and does not take into account the real distributions of mass. Secondly, 
lenses with elliptical deflection potentials were studied (e.g. Kovner 1987b,c; Blandford & Kochanek 
1987a; Kochanek & Blandford 1987; Kochanek et al. 1989). Elliptical potentials are easier to handle than 
elliptical mass distributions. For a small ellipticity elliptical potentials correspond to physical (elliptical) 
mass distributions via the Poisson equation, but for a larger ellipticity (e.g. 0.5 ^ e) the corresponding 
mass distributions obtain unphysical dumbbell shapes (see Kassiola & Kovner 1993). Thirdly, a multipole 
expansion of elliptical mass distributions was considered by Schneider & Weiss (1991) who expanded a 
2-dimensional distribution of mass with elliptical symmetry and calculated the corresponding potential 
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coefScients. Thus, in this approach a (nearly) elliptical mass distribution is described by the several 
lowest-order terms in the expansion. Lastly, for a few special cases the deflection angles were obtained in 
simple (closed) forms. For elliptical “isothermal” lenses (where the surface mass density S ~ [v = 2] 

at large r), the deflection angles were calculated by Kassiola & Kovner (1993) and Kormann, Schneider, & 
Bartelmann (1994). Kassiola & Kovner (1993) also calculated the deflection angles for other integer values 
of V, i.e. = 3, 4, 5, etc. For singular power-law mass distributions, the deflection angle was obtained by 
Grogin & Narayan (1996). 

All of the above lens models avoid numerical integrations and are relatively easy to use. However, 
we know that most (if not all) of them are limited in their applications. For the generalized non-circular 
power-law mass distributions where the parameter v and core radius are arbitrary, numerical integrations 
(Bourassa & Kantowski 1975; Schramm 1990) and multipole expansions (Schneider & Weiss 1991) have 
been considered and the former approach has sometimes been used to calculate lensing properties (e.g. 
most recently by Keeton and Kochanek 1997). In the Schneider & Weiss (1991) approach an elliptical mass 
distribution is expanded and several lowest-order terms are used to describe the original mass distribution. 
This approximation can be effective for a less elliptical mass distribution because only a few terms suffice 
to describe the original mass distribution. However, for a highly elliptical mass distribution this approach 
becomes increasingly less effective because many terms are necessary to approximate the original mass 
distribution. Also, we do not, at present, find any examples in the literature where this approach was used. 
For the above reasons this approach will not be further discussed. In this paper we present an analytical 
approach to lensing by the generalized power-law mass distributions. In our approach all of the terms in 
the infinite series can be integrated in closed forms (as hypergeometric functions), and the series themselves 
have well-defined convergence properties. As a result, we obtain the deflection and magnification in 
ready-to-use forms. Our approach is a straightforward and user-friendly way to calculate lensing properties 
more effectively than the more conventional methods mentioned above. Our results can be used by those 
who want to use (exactly) elliptical surface mass densities (or, 3-dimensional triaxial mass distributions 
directly) and who do not want to use numerical integrations. We applied our results to the Einstein Cross 
and they were very effective (Chae, Turnshek, & Khersonsky 1998). Our approach is as follows. 

We represent the deflection potential as a Fourier series for an arbitrary projected surface density. 
From the deflection potential we obtain the deflection angle as a series. We calculate the projected surface 
mass density of the triaxial power-law mass distribution and thus relate the parameters of the surface 
mass density to those of the 3-dimensional mass density. We calculate the coefficients in the series of the 
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deflection angle for the surface density. We obtain expressions for the deflection angle and magnification 
factor. 

The above procedures are described in §2; the details of mathematics can be found in the Appendices. 
In §3 we describe how the behaviors of critical curves and caustics depend on the parameters of the elliptical 
power-law mass distributions. Concluding remarks are given in §4. 


2. Mathematical Formulation 

2.1. Basic Equations, Definitions and Notation 

In this subsection we briefly summarize the basic equations of gravitational lensing and define the 
coefficients of the series of the deflection potential. Our treatment of lensing is based on the scalar potential 
formalism given by Schneider (1985) and Blandford & Narayan (1986). Our notation is similar to that used 
by Schneider, Ehlers, & Falco (1992, hereafter SEF) and Blandford & Kochanek (1987b). 

When a light ray emitted by a quasar (i.e. the source) at angular diameter distance Ds passes through 
a transparent distribution of mass (i.e. the deflector) at D^, the light ray experiences a deflection and the 
condition for the deflected light ray to reach the observer is given by the lens equation (SEF) 

ri= Dds&iO- ( 1 ) 

Here ff is the position vector of the quasar on the source plane and ^ is the impact vector of the light ray on 
the lens plane, both with respect to the optical axis (defined below). The parameter q;(^) is the deflection 
angle and Dd, Dg and Dds are the angular diameter distances to the deflector, the source, and from the 
deflector to the source, respectively. The optical axis is defined as the infinite line joining the observer and 
the center of mass of the lens. By introducing an arbitrary length scale and defining x = (deflector 
vector), Xs = ff/(^ qDs/D d) (source vector), the above equation becomes the dimensionless lens equation 

Xs = X — d{x), (2) 

where d{x) is the scaled deflection angle which is related to the deflection potential 'i/{x) by d{x) = N't/fx). 
The deflection potential is given in polar coordinates (SEF) by 

-| /‘27r poo 

= — / d(j)' / \n[r'^ + r'"^ — 2rr'cos{(j) — , (3) 

'X Jo Jo 
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where the dimensionless surface density k{x) is obtained by dividing the surface density S(^) by 'Scr, i-e. 
k{x) = T,{^ox)/T,cr, and the critical surface mass density Scr is defined by Ecr = {^'^G/c^)~^{DdDds/Ds)~^■ 

By expanding the logarithmic function in equation (3) for r' > r and r' < r respectively, the deflection 
potential can be written as 


u j.\ , ^cosncj) 
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and the deflection angle components are 
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where the coefficients of the series are defined as follows: 

r27T 


and 


pZTT nr 

^o(^) = / dcf)' / K{r ,(j))r dr ^ 

Jo Jo 
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Jo Jo 

n2'K nr 

Cn{T) = / d(t)' sin n(j)' / n{r\ (f)'dr\ 

Jo Jo 

n2'K nOO ^ 

Dn{r) = J dcj)'cos ncf)' J K{r'dr', 

n2'K nOO ^ 
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n27r nOO 
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From equation (2) the image magnification factor is given by 

I -1 

I I I I 

M = 


det 

dxg 


dx 


( 4 ) 


( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 

(9) 

( 10 ) 
( 11 ) 

( 12 ) 

(13) 

where {^-) is the Jacobian matrix in Xg = Xsix)- A point source at Xg on the source plane forms an image 
at X on the lens plane magnified by a factor of |Af|. 
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2.2. The Projected Surface Mass Density 


The 3-dimensional triaxial power-law distribution of mass is described by 


p{X,Y,Z) = po 


X 


1 + - + ^ + - 


Y 


-vjl 


(14) 


where {X,Y,Z) are the body coordinates attached to the leasing object. The positive constants {a,b,c) 
represent the core sizes along each axis. In passing, we mention that a more conventional parameterization 
is obtained through h = au and c = av, where 1 > u > v by convention. The constant po is the density at 
the center. The radial index ly determines how the mass density decreases at large r. If > 3, the lens has a 
finite total mass (M) and po is related to the total mass by po = M[27ra6cB(3/2, where p,= {v — 3)/2 

and B(a,5) is the Beta function. When = 2, the mass distribution is called “isothermal” since the mass 
density decreases like the singular isothermal sphere at large r. A physical distribution of mass must satisfy 
n > 1 because i/ = 1 corresponds to a constant surface density on the lens plane (see eq. [15]). 


Since a lensing galaxy (or any lensing object) can be oriented in an arbitrary direction relative to the 
lens plane, three parameters are needed to relate the body coordinates {X,Y,Z) to the lens coordinates 
{x,y,z), where z is the direction toward the observer. These three parameters are the Eulerian angles 
{a, (3, 7 ) (see, e.g., Goldstein 1980). We have Xi = ^jaijXj, where Xi = {X,Y,Z), Xj = {x,y,z), and 
the transformation matrix (a^) can be found in Goldstein (1980, eq. [4-47]) with i,j = 1,2,3. Using this 
transformation we find 
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Now the 2-dimensional surface density on the lens plane can be evaluated as 


T,{x,y) 


p{x,y,z)dz 




Po 1 

VA[fix,y)]t^+^' 


( 15 ) 



where f{x,y) = C{x,y) — [B{x,y)]'^/A. Note that /x[= — 3)/2] > —1. 

We rewrite f{x,y) as 

f{x,y) = l + Cxx‘^+Cyy‘^+ Cxyxy, (16) 

where = Ci — BfjA, Cy = C 2 — B^jA, and c^y = 2(C'3 — B 1 B 2 /A). In polar coordinates equation (16) 
becomes 

/(r,(/)) = l + r^[P + (5sin(2(() + S')], (17) 

where P = {c^ + Cy)/ 2 , S = tan“^[(ca: - Cy)/cxy\^ and Q = {c^ - Cy)/{2smS). 


Equations (15) and (17) show that the projected surface density of a triaxial mass distribution (eq. [14]) 
is described by a set of concentric ellipses with constant ellipticity (see, e.g., Stark 1977). The ellipticity of 
an ellipse, e, is given by 


^ min 

e = 1 -= 1 — 


l + |e| 


(18) 


where Tmin and r^ax are the semi-minor and semi-major axes of an ellipse, respectively, and e = QjP■ The 
ratio jej is a measure of the shape of an ellipse and 0 < jej < 1 . 


The dimensionless surface density as a function of the deflector vector is given by 


k (£) = K(r, ^) = 


Ko 


{l + r^[P + Qsm{2(l) + 5')]}'^+!' 
Equation (19a) can be rewritten in the following more convenient form 

\ -(ai-I-1) 


K(r, ^) = Ko ■( 1 -b ( — ) [1 -b e cos 2 ((() - 0 o)] 

. - 


(19a) 


(195) 


where the “core radius” Vc = IIVP and the orientation angle (/)o = f — -f, which becomes the standard 
position angle (P.A., north through east) if e > 0 (see also Keeton & Kochanek 1997). In equations (19a,b) 


22a»+1b(^ + 1,^+1)P0 
\/j4 • Pier 


( 20 ) 


2.3. The Deflection Angle and Magnification Factor 

The coefficients of the deflection angle (equations [7] - [11]) are calculated in Appendix A. Using 
the results of Appendix A, we are now ready to consider the full expressions for the deflection angle and 
magnification factor at x on the lens plane due to the mass distribution given by equation (14) or equation 
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(19). From Appendix A, all odd-numbered coefficients vanish and the even-numbered coefficients can be 
written as follows by introducing three functions of r / 2 m (^)> ^ 2 m(^)]- 


Ao(r) = Ko7rr/(°^(r), 

( ^2m (r) = Ko7rri+2'"4^(r) cos(m(5) 
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( 21 ) 

( 22 ) 


(23) 


Now the deflection angle components, equations (5) and ( 6 ), can be written in the following simple 
forms 
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00 
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Using the results of Appendix A, (r), I^Hir), and can be written as follows. 


For r < Ijy/P - IQI 

00 k 

/(o)(r) = ft.(r)[£2(r)]^^^(-l)'Ci(fc,^,0,/i)[e2(?')]'2Ai(-l - /x; l;£i(r)), 


fe=o 1=0 


for r > 1/ y/P — \Q\ and ^ 7 ^ 0 


00 

/(o)(r) = - ^===-/x(r)[£ 2 (r)]'"^[£ 2 (r)]'= 2 Ai(-fc-/x,-fc-Ai;l;£i(r)) 

r^yp^ — 


k=0 


and for /X = 0 


/(o)(r) = 


r^/P^ - 


: In 


a/P2 _ g2^('p2 _ Q2)^4 2Pr2 -h 1 -f (p2 _ Q^y + p 
V/P2 -Q^ + P 


(24) 

(25) 


(26a) 


(266) 


(26c) 
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(27a) 


and for r > Xj— IQj, if /r is not integer, 
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And for any r > 0 
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In equations (26), (27), and (28) we defined the following: 
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Note that 0 < £i(r) < 1 and 0 < £ 2 (?") < 1 for 0 < r < oo. In equations (26) and (27), the range of r for 


each expression corresponds to the range for which the expression converges rapidly. Outside this range each 


expression will still converge but slowly. Equation (28) converges quickly for any r if to > 2, but converges 
very slowly for to = 1 and {P + IQDr^ << 1. For to = 1 and {P + |(5|)r^ << 1, an alternative form of 
Hmir) can be found in Appendix A (eq. [A9]), which converges very quickly. In cartesian coordinates the 
deflection angle is 

d= aa;i + ay] = {ur cos 4> — a^ sin ^)i + (a^ sin ^ cos (j))]. (34) 



-11 - 


The determinant of the Jacobian matrix (^■) is 
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The magnification factor is then 


M = 


1 

1 ~ (V')ll +V’)22 ) + Ip,11 P’,22 —P’,12 


(45) 


The coefficient functions of the deflection angle and magnification factor, equations (26) - (28) (and 
equations [43] and [44]), contain a hypergeometric function of the same form 2 ^’i(j — k — —k — /r; J + 1; x) 
with k + ^ > —1 and j = 0,1, 2,... The hypergeometric function converges for any \x\ < 1. One can prove 
that the series converges rapidly after about the n-th order, where n « fc + ^ + 1 and higher order terms 
than 2n are almost negligible. Thus, one can evaluate the hypergeometic function easily by truncating the 
series somewhere after 2n, depending on the accuracy desired. Furthermore, once the value of the function 
for a fc is known, the function values for other values of k can be calculated using the Gauss recurrence 
relations (see, e.g., Abramowitz & Stegun 1964). 

For the calculation of light travel time difference (i.e. time delay) between an image pair, say at (ri, (pi) 
and (r 2 , ^ 2 ) on the lens plane, one needs to calculate the difference of F{r) (eq. [12]) at the two points, i.e. 
F{ri) — F{r 2 ) = d(j)' K{r',(j)')\nr' ■ r'dr' as well as the other coefficients (equations [7] - [11]). For 

the expression used to determine the time delay the reader is referred to SEF [equations (5.11) and (5.44)]. 


3. Critical Curves and Caustics of Power-Law Lenses 

Our results can be used to see how the case of an arbitrary value of i' differs from the best-studied 
case of = 2 (i.e. the isothermal lenses studied by Kormann et al. 1994 and Kassiola & Kovner 1993) in 
the structure of critical curves and caustics. The geometry and number of images are determined by the 
structure of the caustics on the source plane. For isothermal lenses, Kormann et al. (1994) showed that five 
different types of caustics are possible in the 2 -dimensional parameter space spanned by axis ratio / (which 
is related to the ellipticity e via e = 1 — /) and core radius. In particular, the simple expressions for the 
deflection enabled them to map out the parameter space in terms of the types of caustics (see their Fig. 9). 
Here we use examples to study the effect of varying radial index v as well as ellipticity and core radius on 
the behavior of critical curves and caustics. 

We compare three different values of v including the isothermal case for the same ellipticities and core 
radii (Figures 1, 2, and 3). This provides us with a qualitative understanding of the difference between 
shallower and steeper distributions of mass. The case of e = 0.7 and Tc = 0.5 in Fig. 1 shows one liplike (i.e. 
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of the diamond shape) caustic which we call the first caustic. For e = 0.5 and Vc = 0.5, the second caustic 
of a liplike (or pseudo-elliptic) shape appears inside the first caustic. For e = 0.3 and Tc = 0.5, two cusps 
of the first caustic are inside the second caustic while the other two are “naked” outside and the second 
caustic is of an elliptical shape. At this stage the second caustic is called “radial” while the first caustic is 
called “tangential” (see, e.g., SEF). A further decrease of e for the same value of rc causes the tangential 
caustic to be completely inside the radial caustic (for e = 0.1 and Vc = 0.5). The last and trivial case is 
that there is no caustic. This happens when we sufficiently increase either e or rc from the case of e = 0.7 
and rc = 0.5 in Fig. 1. In the case of no caustic, multiple imaging is impossible. When there is one liplike 
caustic, three images become possible. When two caustics exist, up to five images are possible. 

For a given value of v the critical behavior can be summarized as follows. For sufficiently large 
ellipticity and core radius, no caustics exist (thereby no critical curves). As we decrease either ellipticity 
or core radius, the first caustic comes into existence. A further decrease gives birth to the second caustic, 
which grows as either ellipticity or core radius decreases, finally enclosing the tangential caustic completely 
for a sufficiently small ellipticity or core radius. 

The effect of varying radial index v on the critical behavior can be qualitatively understood by 
comparing Figures 1, 2, and 3. For example, for e = 0.7 and rc = 0.5, the case u = 1.5 has one liplike 
caustic. For the same ellipticity and core radius, the isothermal case has two caustics. For v = 2.5, two 
cusps of the tangential caustic are inside the radial caustic (i.e. the second caustic is grown relative to the 
first caustic). Now the significance of v on the critical structure is evident. For given ellipticity and core 
radius, a sufficiently shallow profile does not allow any caustic. As the profile steepens, the first caustic will 
appear, then the second caustic, and finally the radial caustic will enclose the tangential caustic completely. 


4. Discussion and Conclusion 

We obtained the complete expressions for the deflection angle and magnification factor as (rapidly) 
converging series for the triaxial mass distribution of equation (14) and for the elliptical surface density of 
equation (19). The calculation of an angle (eq. [34]) at a point on the image plane requires evaluations of 
three functions which are given by equations (26), (27) and (28). For the calculation of the magnification 
(eq. [45]) one needs to evaluate two additional functions (equations [43] and [44]). Since they converge 
very rapidly all these functions can be easily evaluated. The hypergeometric function in equations (26), 
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(27), (28), (43) and (44), which is a well-behaved converging series, can be evaluated rapidly. For instance, 
one can evaluate the hypergeometric function by summing up less than tens of terms. However, since the 
hypergeometric function should be evaluated a number of times (namely, the function should be called 
many times) for the calculation of the deflection and magnification, the execution speed in our approach 
depends on how the hypergeometric function is evaluated. If one evaluates the hypergeometric function 
solely by summing up terms, it typically takes about one milli-second for the calculation of the deflection 
and magnification at a time on a current (SPARC20) unix machine with a fractional uncertainty of about 
10“^ or better. For example, if the deflection angle is one arcsecond, then the calculational error is about 
one hundredth milli-arcsecond or smaller with the execution speed specified above. 

Now that our results have been presented, one can ask how well the mass model (eq. [14] or eq. [19]) 
resembles the true mass distribution of a real lens. An application of our results to the Einstein Cross 
(Q2237-I-0305), which is known to be lensed by an isolated barred spiral galaxy, shows that the model 
is a good approximation to the true mass distribution (Chae et al. 1998). In most of the other lensed 
systems external perturbations are known to be important (Keeton et al. 1997). Thus, unless the external 
perturbations are taken into account properly in fitting observational constraints, it is not possible to 
judge whether or not the mass model is a good description of the true mass distribution. There have been 
numerous lens modeling examples where isothermal models were usedj^ The mass distributions of lenses in 
nature may differ only slightly from isothermal distributions. However, lenses in the universe will not, in 
general, exactly mimic isothermal distributions. Thus, it is worthwhile to consider arbitrary values of v. 
The model we considered here encompasses the isothermal model as a special case (i.e. u = 2) but allows us 
to vary the parameter v from the isothermal value. 

In conclusion, the following are what we consider to be the most notable aspects of our results. 

1. The model we have considered represents the most general case of the familiar power-law mass 
distributions. It can be considered as a generalization of the isothermal model. Hence, our results can 
be used to model various single mass distributions (e.g. elliptical/spiral galaxies) until more physical 
mass models are developed. 

2. The resulting expressions for the deflection angle and magnification factor are exact, and no 
approximations were made in any steps of the calculation. 

^We do not list those examples here because there are many. 
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3. The series expressions for the deflection angle and magnification factor converge rapidly and can be 
easily evaluated with the desired accuracies. 

Appendices 


A. The Series Coefficients of the Deflection Angle 


We evaluate the coefficients of the deflection angle defined in §2.1, i.e., equations (7) - (11), for the 
surface mass distribution of equation (19). All of the odd-numbered coefficients vanish because the surface 
density has symmetry under (j) —>■ (jj + ir. Hence we set n = 2m (m = 1, 2, 3,...). 

We define g{(j)') = P + Qsin{2(j)' + S) for convenience. Equation (8) becomes 


p2-K 

S 2 m(r) = Ko d(j)' cos{2m(j)') 
do 


,2(m+l) 


2{m -|- 1) [1 -I- (/(())')r2]A‘+i 


X 2 E 1 ^ -I- 1,1; m -I- 2; ■ 




(Al) 


1 -I- g{4>')r‘^ 

where 2Ei(a,5;c;x) is the hypergeometric function. We write out the hypergeometric function as a series 
and use the binomial expansion to obtain 
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Jo [l + r‘^P + r‘^Qsm{2(j)'+ S)J+^^+^ 
where (a)t, is the pochhammer symbol, namely (a)t, = r(a -|- 5 )/r(a). The evaluation of the integral in 
equation (A2) can be found in Appendix B. So we have 
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where is the associated Legendre function. The argument 2 : = (1 + r^P)l \/(1 + r'^py — {r^Q)^ 

> 1 and by analytic continuation (Gradshteyn & Ryzhik 1994) it can be shown that 




z-l\ = /z + 1 


lim 
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2Fi{-1 - /i,-l - ^ - m; 7 ; 
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After substituting the above expression into equation (A3), we have 
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Using the defined functions of r (equations [30] and [31]) and the relation — m — = (—l)'"(^ + + +l)m, 


we have 
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As equation (Al) implies, the above expression (eq. [A6]) converges slowly when {P — [QD'c^ becomes 
large, i.e. g{(j)')r‘^/\1 + g((/>0^^] ~ 1) while it converges rapidly for small r. However, equation (Al) can be 
transformed to an alternative form which converges rapidly for r > 1/7 P — jQj- This alternative form is 


given by 
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The coefficient C 2 m(r)/r^™ can be calculated in a very similar way and the resulting expressions are 
identical to equations (A6) and (A7), except that cos[m(7r/2 — S')] is replaced with sin[m(7r/2 — S)]. 


The coefficients D 2 m and E 2 m can also be calculated in a similar way. We find 
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From the first line of equation (A8), it is clear that converges slowly when {P+\Q\)r‘^ « 1. This 

is particularly true for to = 1. For to > 2, the convergence is still acceptable even when (P + |Q|)r^ << 1. 
Thus, we give an alternative form of I^^ir) for to = 1, which converges very rapidly when r < 1 /\/P + jQj: 



r - (-l)fe T{g + k+ 1) , , A /2gy [1 - (-l)^][r(V2 + 1)]^ 
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The above expression can be obtained by dividing the integration interval of r' in equation (10) into [r, ri] 
and [ri, 00] where r < ri < 1/ \/g{ 4 >')- 
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The coefficient ^o(^) is obtained by substituting m=0 into the coefficient B 2 m{'r), i.e. 
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From the hypergeometric function in equation (Al), it is obvious that either for a lower order of m or a 
larger value of r, the series will converge more slowly. Indeed, the coefficient Ao(r) converges very slowly 
for r » Ij^yP — \Q\. However, for large r the coefficient Ao(r) can be represented by another form which 
converges extremely rapidly. This alternative form of Ao(r) can be obtained from equation (7) (for /i 0). 
Equation (7) becomes, after the integration over r'. 
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In the above expression, the second term can be easily calculated using = ^ + ^ + ... for A > 1. So 
we have 
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(All) 


for /X y^ 0 and r > Xj — |g|. When ^ = 0, for any r 
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(A12) 


For a circular distribution of mass on the lens plane (i.e. g = 0) all of the coefficients vanish except for 
Ao{r). 


B. Evaluation of Two Integrals 


For A > |P|(y^ 0) we have 
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In the above equation the second term vanishes. We use the following relationship derived from an integral 
representation of the associated Legendre function (Gradshteyn & Ryzhik 1994) 
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We then have 
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Fig. 1.— Critical Curves and Caustics for the radial index v = 1.5 (see §2.2). Two different core sizes 
(arbitrary scale) and four different ellipticities are considered. For Tc = 0.5, each ellipticity corresponds to a 
different type of caustic, which shows that caustic type evolves as e decreases (see §3). Also note that for a 
fixed ellipticity, caustic type evolves as Tc decreases. 
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Fig. 2.— Critical Curves and Caustics for v = 2 (the “isothermal” distribution). The same core sizes and 
ellipticities as in Fig. 1 are considered. This figure can be compared with Figures 1 and 3 to see how caustic 
type varies as v varies (see §3). 
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Fig. 3.— Critical Curves and Caustics for 1 / = 2.5. The same core sizes and ellipticities as in Fig. 1 are 
considered. This figure can be compared with Figures 1 and 2 to see how caustic type varies as v varies (see 



















